#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>

#define PITEM(p, s, i) ((double *)((uint8_t *)(p) + (s)*(i)))
#define DITEM(p, s, i) (* PITEM(p, s, i))

/** Определение попадания точки в многоугольник.
 *
 * \param[in] x   абсцисса проверяемой точки
 * \param[in] y   ордината проверяемой точки
 * \param[in] n   количество вершин многоугольника
 * \param[in] px  указатель на массив абсцисс многоугольника
 * \param[in] py  указатель на массив ординат многоугольника
 * \param[in] sx  расстояние в байтах между элементами массива абсцисс
 * \param[in] sy  расстояние в байтах между элементами массива ординат
 *
 * \return Функция возвращает 1, если заданная точка попадает в заданный
 *         многоугольник и 0 в противном случае.
 *
 *
 * Из одного форума:
 *
 * \verbatim
 * A> Ребят, никто не знает как можно реализовать попадание точки в
 *    многоугольник за логарифмическое время, но без предобработки?
 * B> т.е. линейное (по числу вершин) не устраивает?
 * \endverbatim
 *
 * Предложенная ниже реализация основывается на материале, найденном по ссылке
 * http://users.livejournal.com/_winnie/237888.html
 */
int p_point_in_polygon(double x, double y,
    int n, const double *px, const double *py, int sx, int sy) {
  int i, parity, below0, below1;
  double x0, y0, x1, y1, det;

  x0 = DITEM(px, sx, n-1) - x; y0 = DITEM(py, sy, n-1) - y;
  below0 = y0 < 0.0;

  for (i = 0, parity = 0; i < n; ++i) {
    x1 = DITEM(px, sx, i) - x; y1 = DITEM(py, sy, i) - y;
    below1 = y1 < 0.0;
    det = (x0*(y1 - y0) - y0*(x1 - x0));
    if (!below0 && below1 && det > 0) ++parity;
    if (below0 && !below1 && det < 0) ++parity;
    x0 = x1; y0 = y1; below0 = below1;
  }

  return parity % 2;
}



/** Рабочий контекст реализации алгоритма Уайлера - Атертона.
 */
typedef struct tag_p_wa_context {
  void *pool;
} p_wa_context_t;


/** Контур для алгоритма Уайлера - Атертона.
 */
typedef struct tag_p_wa_contour {
  int n;
  const double *px;
  const double *py;
  int sx;
  int sy;
  struct tag_p_wa_point *data;
} p_wa_contour_t;


#define VERTEX          0x01
#define VERTEX_INNER    0x03
#define CROSS           0x04
#define CROSS_ENTERING  0x0c



/** Вершина или точка пересечения рёбер для алгоритма Уайлера - Атертона.
 */
typedef struct tag_p_wa_point {
  /** Тип точки. */
  int type;

  double x;
  double y;

  /** Указатель на следующую точку списка. */
  struct tag_p_wa_point *next;

  union {
    struct {  // 0 != (type & VERTEX)
      /** Указатель на следующую ВЕРШИНУ. */
      struct tag_p_wa_point *next_vertex;
    };

    struct {  // 0 != (type & CROSS)
      double t;
      /** Указатель на эту же точку пересечения в другом контуре. */
      struct tag_p_wa_point *cross;
    };
  };
} p_wa_point_t;



/** Добавление точки пересечения к контуру.
 *
 * \param v начальная вершина ребра контура
 * \param c вставляемая в список точка пересечения
 *
 * Использована тривиальная сортировка при вводе; есть основания полагать, что
 * количество элементов будет относительно мало.
 *
 */
static
p_wa_point_t * p_wa_insert_cross(p_wa_point_t *v, p_wa_point_t *c) {
  int parity = (v->type & VERTEX_INNER) == VERTEX_INNER ? 0 : 1;
  while (v && v->next && (v->next->type & CROSS) && (v->next->t < c->t)) {
    v = v->next;
    ++parity;
  }
  c->type = (parity % 2) ? CROSS_ENTERING : CROSS;
  c->next = v->next;
  v->next = c;
  return c;
}


/** Вычисление точки пересечения отрезков.
 */
static
int p_edge_intersect(double x0, double y0, double x1, double y1,
    double X0, double Y0, double X1, double Y1,
    double *cx, double *cy, double *t, double *T) {
  double t1, t2;

  double dx = x1 - x0;
  double dy = y1 - y0;
  double DX = X1 - X0;
  double DY = Y1 - Y0;

  // считаем уравнения прямых проходящих через отрезки
  double a1 = -dy;
  double b1 = +dx;
  double d1 = -(a1*x0 + b1*y0);

  double a2 = -DY;
  double b2 =  DX;
  double d2 = -(a2*X0 + b2*Y0);

  // подставляем концы отрезков, для выяснения в каких полуплоскотях они
  double seg1_line2_start = a2*x0 + b2*y0 + d2;
  double seg1_line2_end =   a2*x1 + b2*y1 + d2;

  double seg2_line1_start = a1*X0 + b1*Y0 + d1;
  double seg2_line1_end =   a1*X1 + b1*Y1 + d1;

  // если концы одного отрезка имеют один знак, значит он в одной полуплоскости
  // и пересечения нет.
  if (seg1_line2_start * seg1_line2_end >= 0
      || seg2_line1_start * seg2_line1_end >= 0) {
    return 0;
  }

  t1 = seg1_line2_start / (seg1_line2_start - seg1_line2_end);
  t2 = seg2_line1_start / (seg2_line1_start - seg2_line1_end);
  *cx = x0 + t1*dx;
  *cy = y0 + t1*dy;
  *t = t1;
  *T = t2;

  return 1;
}


static
p_wa_point_t * p_wa_alloc_vertex(p_wa_context_t *ctx, p_wa_point_t *prev,
    double x, double y, int type) {
  p_wa_point_t *v = malloc(sizeof(p_wa_point_t));  // TODO use pools!
  v->type = type;
  v->x = x; v->y = y;
  v->next = NULL;
  if (prev) prev->next = v;
  v->next_vertex = NULL;
  if (prev) prev->next_vertex = v;
  return v;
}


static
p_wa_point_t * p_wa_alloc_cross(p_wa_context_t *ctx,
    double x, double y, double t) {
  p_wa_point_t *v = malloc(sizeof(p_wa_point_t));  // TODO use pools!
  v->x = x; v->y = y;
  v->next = NULL;
  v->t = t;
  return v;
}


static
int p_wa_point_in_polygons(double x, double y, int N, p_wa_contour_t *cc) {
  int i, parity;
  p_wa_contour_t *c;
  for (parity = 0, i = 0; i < N; ++i) {
    c = cc + i;
    parity += p_point_in_polygon(x, y, c->n, c->px, c->py, c->sx, c->sy);
  }
  return parity % 2;
}


static
p_wa_point_t * p_wa_alloc_contour_vertices(p_wa_context_t *ctx,
    p_wa_contour_t *contour, int N, p_wa_contour_t *another_poly_contours) {
  int i;
  double x, y;
  p_wa_point_t *v0, *v;

  for (v = NULL, i = 0; i < contour->n; ++i) {
    x = DITEM(contour->px, contour->sx, i);
    y = DITEM(contour->py, contour->sy, i);
    v = p_wa_alloc_vertex(ctx, v, x, y,
        p_wa_point_in_polygons(x, y, N, another_poly_contours)?
        VERTEX_INNER : VERTEX);
    if (i == 0) v0 = v;
  }
  if (v != NULL) {
    v->next = v0;
    v->next_vertex = v0;
  }

  return v0;
}


static
void p_wa_print_contour(const p_wa_point_t *start, int level, int depth) {
  const char *type;
  int i;
  double x1, y1, x2, y2;
  const p_wa_point_t *v = start;

  if (level > depth) return;

  do {
    if ((v->type & VERTEX_INNER) == VERTEX_INNER) {
      type = "inner vertex";
    } else if ((v->type & VERTEX) == VERTEX) {
      type = "outer vertex";
    } else if ((v->type & CROSS_ENTERING) == CROSS_ENTERING) {
      type = "entering cross";
    } else if ((v->type & CROSS) == CROSS) {
      type = "exiting cross";
    } else {
      type = "UNKNOWN";
    }

    for (i = 0; i < level; ++i) printf("  ");
    printf("(%g, %g): %s", v->x, v->y, type);

    if (v->type & CROSS) {
      if (v->cross == NULL) {
        printf(" -> NULL\n");
      } else {
        x1 = v->next->x - v->x;
        x2 = v->cross->next->x - v->x;
        y1 = v->next->y - v->y;
        y2 = v->cross->next->y - v->y;
        printf(" -> %s\n", x1*y2 - x2*y1 > 0.0 ? "left" : "right");
        p_wa_print_contour(v->cross->next, level + 1, depth);
      }
    } else {
      printf("\n");
    }

    v = v->next;
  } while (v && v != start);
}


static
void p_wa_trace_contour(const p_wa_point_t *start) {
  double x1, y1, x2, y2;
  int i = 0;

  const p_wa_point_t *v = start;
  do {
    v = v->next;
  //} while (v && (v->type & CROSS) != CROSS && v != start);
  } while (v && (v->type & CROSS_ENTERING) != CROSS_ENTERING && v != start);

  start = v;
  do {
    printf("%04d: (%g, %g)\n", ++i, v->x, v->y);

    if (v->type & CROSS) {
      x1 = v->next->x - v->x;
      x2 = v->cross->next->x - v->x;
      y1 = v->next->y - v->y;
      y2 = v->cross->next->y - v->y;
      if (x1*y2 - x2*y1 < 0) {
        v = v->cross->next;
      } else {
        v = v->next;
      }
    } else {
      v = v->next;
    }
  } while (v != start && v->cross != start);  // можем подойти с разных сторон!
}


/** */
void p_wa_intersect(int N1, p_wa_contour_t *P1, int N2, p_wa_contour_t *P2) {
  p_wa_context_t ctx;
  int i, j;
  p_wa_contour_t *c1, *c2;
  p_wa_point_t *p1, *p2, *C1, *C2;
  double cx, cy, t1, t2;

  // создаём списки вершин для каждого контура первого полигона
  for (i = 0; i < N1; ++i) {
    P1[i].data = p_wa_alloc_contour_vertices(&ctx, P1+i, N2, P2);
  }

  // создаём списки вершин для каждого контура второго полигона
  for (i = 0; i < N2; ++i) {
    P2[i].data = p_wa_alloc_contour_vertices(&ctx, P2+i, N1, P1);
  }

  // ищем попарные пересечения контуров
  for (i = 0; i < N1; ++i) {
    c1 = P1 + i;
    for (j = 0; j < N2; ++j) {
      c2 = P2 + j;

      p1 = c1->data;
      do {
        p2 = c2->data;
        do {
          if (p_edge_intersect(
                p1->x, p1->y, p1->next_vertex->x, p1->next_vertex->y,
                p2->x, p2->y, p2->next_vertex->x, p2->next_vertex->y,
                &cx, &cy, &t1, &t2)) {
            C1 = p_wa_insert_cross(p1, p_wa_alloc_cross(&ctx, cx, cy, t1));
            C2 = p_wa_insert_cross(p2, p_wa_alloc_cross(&ctx, cx, cy, t2));
            C1->cross = C2; C2->cross = C1;
          }
        } while ((p2 = p2->next_vertex) != c2->data);
      } while ((p1 = p1->next_vertex) != c1->data);
    }
  }

  puts("*** POLYGON P1 ***");
  for (i = 0; i < N1; ++i) {
    p_wa_print_contour(P1[i].data, 0, 1);
  }

  puts("*** POLYGON P2 ***");
  for (i = 0; i < N2; ++i) {
    p_wa_print_contour(P2[i].data, 0, 1);
  }

  p_wa_trace_contour(P1[0].data);
}


int main() {
  double x1[] = { -1,  1,  1, -1 };
  double y1[] = {  2,  2, -2, -2 };
  double x2[] = {  0,  2,  2,  0 };
  double y2[] = {  1,  1, -1, -1 };

  p_wa_contour_t C1 = {
    sizeof(x1)/sizeof(x1[0]), x1, y1, sizeof(x1[0]), sizeof(y1[0]), NULL
  };

  p_wa_contour_t C2 = {
    sizeof(x2)/sizeof(x2[0]), x2, y2, sizeof(x2[0]), sizeof(y2[0]), NULL
  };

  p_wa_intersect(1, &C1, 1, &C2);

  return 0;
}

